Вы узнаете

  • Как можно смоделировать "поведение" дисперсии

Вы сможете

  • Построить модели для данных, демонстрирующих высокую степень гетерогености дисперсии
  • Подобрать оптимальную модель, учитывающую ковариаты дисперсии
  • Построить смешанную модель учитывающую не только группирующие (случайные) факторы, но и гетерогенность дисперсии.

В начале курса мы записывали "обычную" регрессионную модель в таком виде

\[y_i = \beta_0 + \beta_1x_i + \epsilon_i\]

Фиксированная часть модели: \(y_i = \beta_0 + \beta_1x_i\)
Случайная часть модели: \(\epsilon_i\)

Важное ограничивающее условие метода: \(\epsilon \sim N(0, \sigma^2)\)

Теория смешанных моделей расширяет наши представления о случайной части модели

из Zuur et al. 2009

Модели на языке матриц

Простая линейная модель \[y_i = \mathbf X_i\mathbf\beta + \epsilon_i\]

\[\epsilon_i \sim N(0, \mathbf\sigma^2_i)\]

Смешанная линейная модель с группирующими факторами

\[y_i = \mathbf X_i\mathbf\beta + \mathbf Z_i\mathbf b_i + \epsilon_i\] \[\epsilon \sim N(0, \Sigma_i)\] \[b_i \sim N(0, \mathbf\Psi)\]

Расширенная смешанная линейная модель \[y_i = \mathbf X_i\mathbf\beta + \mathbf Z_i\mathbf b_i + \epsilon_i\] \[\epsilon \sim N(0, \sigma^2 \times \mathbf\Lambda_i)\] \[b_i \sim N(0, \mathbf\Psi)\]

Моделирование гетерогенности дисперсии

Способствуют ли взрослые мидии притоку молоди?

Даные взяты из работы Khaitov, 2013

myt <- read.table("data/myt.csv", sep = ";", header = T)
head(myt, 12)
##    Year Bank Sample Recruits Small Large
## 1  1997 vor2      1        0    12    25
## 2  1997 vor2      2        0    23    27
## 3  1997 vor2      3        0    17    36
## 4  1997 vor2      4        0    17    27
## 5  1997 vor2      5        0    30    25
## 6  1997 vor2      6        0    12    27
## 7  1999 vor2      1       12    12    41
## 8  1999 vor2      2        2    16    38
## 9  1999 vor2      3        1    13    49
## 10 1999 vor2      4        6    19    52
## 11 1999 vor2      5        1     7    31
## 12 1999 vor2      6        2     9    34

В качестве зависимой перменной будем анализировать \(\sqrt{N_{recruits}}\)

myt$Sq_Recruits <- sqrt(myt$Recruits)

myt$fYear <- factor(myt$Year)

Строим обычную регрессионную модель

mod_formula <- formula(Sq_Recruits ~ Large + fYear + Bank + Large:fYear + 
    Large:Bank)

M1_lm <- lm(mod_formula, data = myt)

anova(M1_lm)
## Analysis of Variance Table
## 
## Response: Sq_Recruits
##              Df Sum Sq Mean Sq F value  Pr(>F)    
## Large         1    481     481   87.41 < 2e-16 ***
## fYear        13   2189     168   30.59 < 2e-16 ***
## Bank          2    344     172   31.27 1.2e-12 ***
## Large:fYear  13    228      18    3.18 0.00021 ***
## Large:Bank    2      2       1    0.14 0.86594    
## Residuals   217   1194       6                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Можем ли мы доверять этим результатам?

Проведем диагностику данной модели

Мы не можем доверять величинам уровня значимости! Нарушается условие применимости линейных моделей, основанных на нормальном распределении остатков с неизменной дисперсией.

Ковариата дисперсии (Variance covariate)

Видно, что дисперсия остатков убывает по мере увеличения значения переменной Large и меняется год от года

Переменные Large и fYear являются ковариатами дисперсии

Ковариата дисперсии (Variance covariate)

Новая модель будет включен еще один компонент \[\epsilon \sim N(0, \sigma^2 \times f(VC))\]

\(VC\) - ковариата дисперсии
\(f(VC)\) - функция, вводящая попаравку, стабилизирующую дисперсию
В зависимости от формы функции \(f(VC)\) мы получим разную структуру дисперсии в модели

Различные формы структуры дисперсии

Для дальнейших вычислений необходимо использовать функцию gls из пакета nlme

library(nlme)
M1_gls <- gls(mod_formula, data = myt)

Эта функция дает результаты полностью идентичные результатам функции lm()

anova(M1_gls)
## Denom. DF: 217 
##             numDF F-value p-value
## (Intercept)     1    1408  <.0001
## Large           1      87  <.0001
## fYear          13      31  <.0001
## Bank            2      31  <.0001
## Large:fYear    13       3  0.0002
## Large:Bank      2       0  0.8659

Фиксированная структура дисперсии

Предполагается, что дисперсия изменяется пропорционально значениям ковариаты дисперсии

\[\epsilon_i \sim N(0, \sigma^2 \times VC_i)\]

В данном случае нет небходимости в подборе неизвестных параметров для функции \(f(VC)\)

M2_gls <- gls(mod_formula, data = myt, weights = varFixed(~Large))

Сравним две модели

AIC(M1_gls, M2_gls)
##        df  AIC
## M1_gls 33 1228
## M2_gls 33 1290

НО! Дисперсия явно убывала по мере увеличения значений перменной Large

Разные дисперсии для разных уровней категориальных предикторов

\[\epsilon_{i,j} \sim N(0, \sigma^2_j)\]

При построении моделей с такой структурой дисперсии подбирается M-1 новых параметров, где M - количество уровней категориального предиктора.

M3_gls <- gls(mod_formula, data = myt, weights = varIdent(form = ~1 | fYear))

Сравнение моделей

Важно! Модели M1_gls и M3_gls вложенные

M1_gls: \(\sigma^2_1 = \sigma^2_2 = ... = \sigma^2_m\)
M3_gls: \(k_1\sigma^2_1 = k_2\sigma^2_2 = ... = k_m\sigma^2_m\)

anova(M1_gls, M3_gls)
##        Model df  AIC  BIC logLik   Test L.Ratio p-value
## M1_gls     1 33 1228 1340   -581                       
## M3_gls     2 46 1131 1287   -520 1 vs 2     123  <.0001

Модель M3_gls лучше!

У нас два категориальных предиктора

M3_gls2 <- gls(mod_formula, data = myt, weights = varIdent(form = ~1 | 
    Bank))

anova(M1_gls, M3_gls2)
##         Model df  AIC  BIC logLik   Test L.Ratio p-value
## M1_gls      1 33 1228 1340   -581                       
## M3_gls2     2 35 1216 1334   -573 1 vs 2    16.5  0.0003

Что произошло в результате работы функции varIdent()?

summary(M3_gls)

Часть вывода summary(M3_gls)

Variance function:    
 Structure: Different standard deviations per stratum     
 Formula: ~1 | fYear     
 Parameter estimates:     
 1997  1999  2000  2001  2002  2003  2004  2005  2006  2007  2008      
 1.00  2.62  4.39  3.47  2.84  5.85  4.93  3.21  2.95  3.87  7.98      
 2009  2010  2011     
 9.26  5.97 13.59

Степенная зависимость дисперсии от ковариаты

\[\epsilon_{i,j} \sim N(0, \sigma^2 \times |VC|^{2\delta})\]

Параметр \(\delta\) неизвестен и требует оценки

Если \(\delta = 0\), то струкутра дисперсии будет аналогична структуре диcперсии в "обычной" регрессионной модели, где \(\epsilon \sim N(0, \sigma^2)\)

Важно! Если значения ковариаты дисперсии могут принимать значение равное нулю, то такая форма струкутры дисперсии неопределена и использоваться не может.

M4_gls <- gls(mod_formula, data = myt, weights = varPower(form = ~Large))

Степенная зависимость дисперсии от ковариаты

Оценка параметра \(\delta\)

M4_gls$modelStruct
## varStruct  parameters:
##  power 
## -0.214

Задание

Степенная зависимость дисперсии от ковариаты может учитывать и взаимодействие ковариаты дисперсии с категориальными предикторами

### Напишите код, с помощью которого в модели будет учтена степенная зависимость дисперсии от переменной Large, но разная для каждого уровня фактора fYear. Аналогичный код напишите для фактора Bank

Hint Изучите справку по функции varPower()

Решение

M5_gls <- gls(mod_formula, data = myt, weights = varPower(form = ~Large | 
    fYear))
M6_gls <- gls(mod_formula, data = myt, weights = varPower(form = ~Large | 
    Bank))
M5_gls$modelStruct
## varStruct  parameters:
##    1997    1999    2000    2001    2002    2003    2004    2005 
## -0.4236 -0.3033 -0.1106 -0.1966 -0.2512 -0.0680 -0.1480 -0.2494 
##    2006    2007    2008    2009    2010    2011 
## -0.2805 -0.1650  0.0359  0.0860 -0.0536  0.1860
M6_gls$modelStruct
## varStruct  parameters:
##    vor2    vor4    vor5 
## -0.1980 -0.0488 -0.2299

Экспоненциальная зависимость дисперсии от ковариаты

\[\epsilon_{i,j} \sim N(0, \sigma^2 \times e^{2\delta \times VC_i})\]

Эта форма структуры дисперсии может применяться для случаев, когда \(VC = 0\)

Если \(\delta = 0\), то структура дисперсии будет аналогична струкуте диспесии в "обычной" регрессионной модели, то есть \(\epsilon_{i,j} \sim N(0, \sigma^2)\)

M7_gls <- gls(mod_formula, data = myt, weights = varExp(form = ~Large))
M8_gls <- gls(mod_formula, data = myt, weights = varExp(form = ~Large | 
    fYear))
M9_gls <- gls(mod_formula, data = myt, weights = varExp(form = ~Large | 
    Bank))

Экспоненциальная зависимость дисперсии от ковариаты

Оцененные параметры

M7_gls$modelStruct
## varStruct  parameters:
## expon 
## -0.02
M8_gls$modelStruct
## varStruct  parameters:
##     1997     1999     2000     2001     2002     2003     2004 
## -0.03093 -0.02296 -0.00946 -0.01591 -0.01679 -0.00395 -0.01906 
##     2005     2006     2007     2008     2009     2010     2011 
## -0.02190 -0.02417 -0.01310  0.00743  0.00939 -0.00347  0.02475
M9_gls$modelStruct
## varStruct  parameters:
##     vor2     vor4     vor5 
## -0.02360 -0.00896 -0.02477

Усложненная степенная зависимость дисперсии от ковариаты

\[\epsilon_{i,j} \sim N(0, \sigma^2 \times (\delta_1 + |VC|^{2\delta_2})^2)\]

Впорос:

При каких значениях параметров функции \(f(VC)\) cтруктура дисперсии будет аналогична структуре дисперсии в "обычной" регрессионной модели?

Ответ

При \(\delta_1=0\) и \(\delta_2=0\) выражение \(\epsilon_{i,j} \sim N(0,\sigma^2 \times (0 + |VC|^{0})\) будет эквивалентно \(\epsilon_{i,j} \sim N(0, \sigma^2)\)

Усложненная степенная зависимость дисперсии от ковариаты

M10_gls <- gls(mod_formula, data = myt, weights = varConstPower(form = ~Large))
# M11_gls <-gls(mod_formula, data = myt, weights = varConstPower(form =
# ~ Large|fYear))
M12_gls <- gls(mod_formula, data = myt, weights = varConstPower(form = ~Large | 
    Bank))

Оцененные параметры

M10_gls$modelStruct
## varStruct  parameters:
##   const   power 
## -17.023  -0.214
M12_gls$modelStruct
## varStruct  parameters:
## const.vor2 const.vor4 const.vor5 power.vor2 power.vor4 power.vor5 
##    -1.3195   -17.2797    -1.9165    -0.2939     0.0127    -0.2608

Комбинированная структура дисперсии

M13_gls <- gls(mod_formula, data = myt, weights = varComb(varIdent(form = ~fYear), 
    varPower(form = ~Large)))
M14_gls <- gls(mod_formula, data = myt, weights = varComb(varIdent(form = ~Bank), 
    varPower(form = ~Large)))
M15_gls <- gls(mod_formula, data = myt, weights = varComb(varIdent(form = ~fYear), 
    varExp(form = ~Large)))
M16_gls <- gls(mod_formula, data = myt, weights = varComb(varIdent(form = ~Bank), 
    varExp(form = ~Large)))

Задание

Найдите модель с наилучшей структурой дисперсии

Решение

AICs <- AIC(M1_gls, M2_gls, M3_gls, M4_gls, M5_gls, M6_gls, M7_gls, M8_gls, 
    M9_gls, M10_gls, M12_gls, M13_gls, M14_gls, M15_gls, M16_gls)

Решение

AICs[AICs$AIC == min(AICs$AIC), ]
##        df  AIC
## M5_gls 47 1131
M5_gls$call
## gls(model = mod_formula, data = myt, weights = varPower(form = ~Large | 
##     fYear))

Диагностика модели с оптимальной структурой дисперсии

Диагностика модели с оптимальной структурой дисперсии

Диагностика модели с оптимальной структурой дисперсии

Диагностика модели с оптимальной структурой дисперсии

Можно ли упростить модель?

M5_gls_ML <- update(M5_gls, method = "ML")
# С какого-то момента перестал работать drop1() drop1(M5_gls_ML, test =
# 'Chi')
M5_gls_ML_a <- update(M5_gls_ML, . ~ . - Large:fYear)
M5_gls_ML_b <- update(M5_gls_ML, . ~ . - Large:Bank)
anova(M5_gls_ML, M5_gls_ML_a)
##             Model df  AIC  BIC logLik   Test L.Ratio p-value
## M5_gls_ML       1 47 1044 1210   -475                       
## M5_gls_ML_a     2 34 1048 1168   -490 1 vs 2      30  0.0047
anova(M5_gls_ML, M5_gls_ML_b)
##             Model df  AIC  BIC logLik   Test L.Ratio p-value
## M5_gls_ML       1 47 1044 1210   -475                       
## M5_gls_ML_b     2 45 1050 1208   -480 1 vs 2    9.68  0.0079

Эту модель упростить нельзя!
То есть, изменение структуры дисперсии заставляет формулировать иные биологические выводы.

Структура дисперсии может иметь определенный биологический смысл

qplot(x = c(1997, 1999:2011), y = as.vector(unlist(M5_gls$modelStruct))) + 
    xlab("Годы") + ylab("Delta")

  • В большинстве случаев параметр \(\delta\) < 0
  • чем больше обилие взрослых мидий, тем меньше варьирует обилие молоди
  • Есть какая-то многолетняя динамика влияния обилия взрослых на "пятнистость" распределения молоди

Моделирование структуры дисперсии при наличии группирующих (случайных) факторов

Рост крыс при разной диете

Пример взят из книги Pinheiro & Bates, 2000 (Hand and Crowder (1996))

Три группы крыс, содержались при разных условиях кормления 64 дня. Каждую крысу взвешивали с определнной периодичностью.

data("BodyWeight")

bw <- as.data.frame(BodyWeight)

head(bw, 14)
##    weight Time Rat Diet
## 1     240    1   1    1
## 2     250    8   1    1
## 3     255   15   1    1
## 4     260   22   1    1
## 5     262   29   1    1
## 6     258   36   1    1
## 7     266   43   1    1
## 8     266   44   1    1
## 9     265   50   1    1
## 10    272   57   1    1
## 11    278   64   1    1
## 12    225    1   2    1
## 13    230    8   2    1
## 14    230   15   2    1

Задание:

Постройте модель, которая дала бы ответ на вопрос: Изменяется ли характер роста крыс в зависимости от типа диеты?

Решение

M1 <- gls(weight ~ Time * Diet, data = bw)

Все ли хорошо?

  • Важно! Строить простую линейную модель в данном случае некорректо!
  • Дизайн эксперимента изначально включает случайный фактор. Здесь мы имеем дело с повторными наблюдениями одного и того же объекта.
  • Однако мы рассмотрим M1 для демонстрации того, что происходит если не учитывать этой осбенности экспериментального дизайна.

Ризультаты для неправильной модели

anova(M1)
## Denom. DF: 170 
##             numDF F-value p-value
## (Intercept)     1   22268  <.0001
## Time            1      20  <.0001
## Diet            2    1114  <.0001
## Time:Diet       2       2    0.17

Случайные факторы

Важно! в этом эксперименте присутсвует случайный (группирующий) фактор Rat, который необходимо учесть в модели

M2 <- lme(weight ~ Time * Diet, data = bw, random = ~1 | Rat)
M3 <- lme(weight ~ Time * Diet, data = bw, random = ~1 + Time | Rat)

Какую из моделей выбрать?

AIC(M1, M2, M3)
##    df  AIC
## M1  7 1739
## M2  8 1248
## M3 10 1172

Пытаемся ответить на вопрос

anova(M3)
##             numDF denDF F-value p-value
## (Intercept)     1   157    1713  <.0001
## Time            1   157      83  <.0001
## Diet            2    13      85  <.0001
## Time:Diet       2   157       8  0.0007

Рассеяние остатков в модели

Моделируем структуру дисперсии

M3_1 <- update(M3, weights = varIdent(form = ~1 | Diet))

M3_2 <- update(M3, weights = varPower(form = ~Time))

M3_3 <- update(M3, weights = varPower(form = ~Time | Diet))

M3_4 <- update(M3, weights = varExp(form = ~Time))

M3_5 <- update(M3, weights = varExp(form = ~Time | Diet))

M3_6 <- update(M3, weights = varComb(varExp(form = ~Time), varIdent(form = ~1 | 
    Diet)))

Выбираем лучшую модель

AIC(M3, M3_1, M3_2, M3_3, M3_4, M3_5, M3_6)
##      df  AIC
## M3   10 1172
## M3_1 12 1164
## M3_2 11 1173
## M3_3 13 1158
## M3_4 11 1174
## M3_5 13 1155
## M3_6 13 1165

Отвечаем на вопрос

anova(M3_5)
##             numDF denDF F-value p-value
## (Intercept)     1   157    1705  <.0001
## Time            1   157      83  <.0001
## Diet            2    13      85  <.0001
## Time:Diet       2   157       9  0.0003

Смотрим на предсказания модели

MyData <- expand.grid(Time = 1:64, Diet = factor(1:3))


MyData$Predicted <- predict(M3_5, newdata = MyData, level = 0)

ggplot(MyData, aes(x = Time, y = Predicted, color = Diet)) + geom_line(size = 1.5) + 
    geom_point(data = bw, aes(x = Time, y = weight), position = position_jitter())

Summary

Проблему гетерогенности дисперсии можно решить двумя способами:
1.Преобразование перменных
2.Введение в модель той или иной структуры дисперсии, учитывающей тот или иной набор ковариат дисперсии.

Что почитать

  • Zuur, A.F. et al. 2009. Mixed effects models and extensions in ecology with R. - Statistics for biology and health. Springer, New York, NY.

  • Pinheiro J, Bates D (2000) Mixed effects models in S and S-Plus. Springer-Verlag, New York, USA